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Abstract 

Using Hermite's formulation of polynomial stability conditions, static output 
feedback (SOF) controller design can be formulated as a polynomial matrix in- 
equality (PMI), a (generally nonconvex) nonlinear semidefinite programming prob- 
lem that can be solved (locally) with PENNON, an implementation of a penalty 
method. Typically, Hermite SOF PMI problems are badly scaled and experiments 
reveal that this has a negative impact on the overall performance of the solver. 
In this note we recall the algebraic interpretation of Hermite's quadratic form as 
a particular Bezoutian and we use results on polynomial interpolation to express 
the Hermite PMI in a Lagrange polynomial basis, as an alternative to the conven- 
tional power basis. Numerical experiments on benchmark problem instances show 
the substantial improvement brought by the approach, in terms of problem scaling, 
number of iterations and convergence behavior of PENNON. 

Keywords: Static output feedback, Hermite stability criterion, Polynomial matrix inequality, 

Nonlinear semidefinite programming. 



1 Introduction 

In 1854 the French mathematician Charles Hermite studied quadratic forms for counting 
the number of roots of a polynomial in the upper half of the complex plane (or, by a 
simple rotation, in the left half-plane), more than two decades before Routh, who was 
apparently not aware of Hermite's work, see [8]. Hurwitz himself used some of Hermite's 
ideas to derive in 1895 his celebrated algebraic criterion for polynomial stability, now 
called the Routh-Hurwitz criterion and taught to engineering students in tabular form. 
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Hermite's criterion can be interpreted as a symmetric formulation of the Routh-Hurwitz 
criterion. This symmetry can be exploited in a semidefmite programming framework, as 
shown in [3] and [J| in the context of simultaneous stabilization of linear systems. Along 
the same vein, in [5] the problem of static output feedback (SOF) design was formulated 
as a polynomial matrix inequality (PMI) problem. In some cases (e.g. only one input or 
output available for feedback) this PMI problem simplifies to a bilinear matrix inequality 
(BMI) that can be solved numerically with PENBMI, a particular instance of PENNON, 
a general penalty method for nonlinear and semidefinite programming. Only convergence 
to a local optimum is guaranteed, but experiments reported in [5] show that quite often 
the approach is viable numerically. In particular, the SOF PMI formulation involves only 
controller parameters, and does not introduce (a typically large number of) Lyapunov 
variables. 

Our motivation in this paper is to contribute along the lines initiated in [5] and to study 
the impact of SOF PMI problem formulation on the behavior of PENNON, in particular 
w.r.t. data scaling and number of iterations. The Hermite matrix depends quadratically 
on coefficients of the characteristic polynomial, in turn depending polynomially on the 
controller parameters. As a result, coefficients of a given Hermite matrix typically dif- 
fer by several orders of magnitude, and experiments reveal that this poor data scaling 
significantly impacts on the performance of PENNON. 

In this paper we use an alternative formulation of the Hermite matrix, using a Lagrange 
polynomial basis instead of the standard power basis. We build on previous work from the 
computer algebra and real algebraic geometry communities, recalling the interpretation 
of Hermite's quadratic form as a particular Bezoutian, the resultant of two polynomials, 
see [6] and references therein. This interpretation provides a natural choice for the nodes 
of the Lagrange basis. The construction of the Hermite matrix in this basis is carried out 
efficiently by interpolation, overcoming difficulties inherent to Vandermonde matrices, as 
suggested in [12] for general Bezout matrices. 

In addition to digesting and tailoring to our needs results from computational algebraic 
geometry, another contribution of our paper is to extend slightly the characterization of 
[T2] to Hermitian forms with complex and repeated interpolation nodes. In particular, 
in our SOF design application framework, these nodes are roots of either imaginary or 
real part of a target characteristic polynomial featuring spectral properties desirable for 
the closed-loop system. This target polynomial is the main tuning parameter of our 
approach, and we provide numerical evidence that a suitably choice of target polynomial, 
compatible with achievable closed-loop dynamics, results in a significant improvement 
of SOF PMI problem scaling, with positive effects on the overall behavior (convergence, 
number of outer and inner iterations, linesearch steps) of PENNON. Furthermore, some 
of the problems that were not solvable in the power basis, see [5], can now be solved in 
the Lagrange basis. These improvements are illustrated on numerical examples extracted 
from the publicly available benchmark collection COMPLjb, see [9]. 
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2 PMI formulation of SOF design problem 



We briefly recall the polynomial matrix inequality (PMI) formulation of static output 
feedback (SOF) design problem proposed in [5]. 

Consider the linear system 

x = Ax + Bu 
y = Cx 

of order n with m inputs and p outputs, that we want to stabilize by static output feedback 
(SOF) 

u = Ky. 

In other words, given matrices A £ M nxn , B £ M nxm , C £ R pxn , we want to find matrix 
K £ M mxp such that the eigenvalues of closed-loop matrix A + BKC all belong to the 
left half of the complex plane. 

Let k £ M. mp be the vector obtained by stacking the columns of matrix K. Define 

n 

q(s, k) = det (si — A — BKC) = ^ q i (k)s i (1) 

i=0 

as the characteristic polynomial of matrix A + BKC. Coefficients of increasing powers of 
indeterminate s in polynomial q(s, k) are multivariate polynomials in k, i.e. 

ft(*) = X)««*° ( 2 ) 

a 

where a £ N mp describes all monomial powers. 

The Routh-Hurwitz criterion for stability of polynomials has a symmetric version called 
the Hermite criterion. A polynomial is stable if and only if its Hermite matrix, quadratic 
in the polynomial coefficients, is positive definite. Algebraically, the Hermite matrix can 
be defined via the Bezoutian, a symmetric form of the resultant. 

Let a(u), b(u) be two polynomials of degree n of the indeterminate u. Define the bivariate 
quadratic form 

a(u)b(v) - a(v)b(u) i-i j-i 

i=i j=i 

The n-by-n matrix with entries bij is the Bezoutian matrix, whose determinant is the 
resultant of a and b, obtained by eliminating variable u from the system of equations 
a(u) = b(u) = 0. 

The Hermite matrix in power basis of q(s,k), denoted by H p (k), is defined as the 
Bezoutian matrix of the real and imaginary parts of q(ju, k): 

a(u,k) = lmq(jw,k) 
b(u,k) = Req(jw,k). 

The roots of polynomial q(s, k) belongs to the left half-plane if and only if 

n n 
i=0 3=0 



The above relation is a matrix inequality depending polynomially on parameters k. There- 
fore, finding k amounts to solving a polynomial matrix inequality (PMI) problem. 

Example 2.1 As an illustrative example, consider problem NN6 in The closed-loop 
characteristic polynomial is (to 8 significant digits): 

q{s, k) = s 9 + 23.300000s 8 + (4007.6500 - 14.688300fc 2 + 14.685000fc 4 )s 7 
+ (91133.935 - 14.685000A;i + 14.688300fc 3 + 15.132810£;4)s 6 
+ (1149834.9 - 57334.489£; 2 + 15.132810^ + 36171.693A;4)s 5 
+ (20216420 - 57334.489^ + 36171.693*3 + 35714.763fc 4 )s 4 
+ (49276365 - 12660338A; 2 + 35714.763Jfc 3 + 3174671. 8A; 4 )s 3 
+ (-1562.6281 • 10 5 - 12660338*1 - 3174671. 8k 3 + 3133948.9A;4)s 2 
+ (-4315.5562 • 10 5 + 95113415£; 2 + 3133948.9Jfe 3 )s 
+ 95113415*1 

with SOF gain K = [k% k 2 k 3 k^] . The 9-by-9 Hermite matrix of this polynomial cannot 
be displayed entirely for space reasons, so we choose two representative entries: 

if£ 3 (fc) = 10244466 ■ 10 8 - 53923375 ■ 10 7 h + 55487273 ■ 10% 

+ 10310826 • 10 7 /fc 3 - 32624061 • 10 7 * 4 + 16028416 • lO 7 ^^ 

- 27103829 • I0 4 hk 3 - 36752006 ■ 10%fc 4 

- 43632833 ■ lO 6 *;^ - 43073807 ■ lO^hh 

+ 22414163/^ + 10078541 • lO 6 **^ + 99492593 • 10 5 ^ 

and 

H[ g (k) = 23.300000. 

We observe that this Hermite matrix is ill-scaled, in the sense that the coefficients of 
its entries (multivariate polynomials in kj) differ by several orders of magnitude. This 
representation is not suitable for a matrix inequality solver. 

3 A simple scaling strategy 

A possible remedy to address the poor scaling properties of the Hermite matrix is to scale 
the frequency variable s, that is, to substitute ps for s in the characteristic polynomial 
q(s,k), for a suitable positive scaling p. Finding the optimal value of p (e.g. in terms 
of relative scaling of the coefficients of the Hermite matrix) may be formulated as an 
optimization problem, but numerical experiments indicate that nearly optimal results are 
achieved when following the basic strategy consisting of choosing p such that the constant 
and highest power polynomial coefficients are both equal to one. For example, this idea 
was implemented by Huibert Kwakernaak in the scale function of the Polynomial Toolbox 
for Mat lab, see [TT] . 
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Example 3.1 Consider the simple example AC4 in [9j. The open-loop characteristic 
polynomial is 

q{s,0) = det(sl - A) 

= s 4 + 150.92600s 3 + 130.03210s 2 - 1330.6306s - 66.837750 

with Hermite matrix in power basis 

" 88936.354 10087.554 

-162937.14 1330.631 

10087.554 20955.855 

1330.6306 150.92600 



H 1 



To measure quantitatively the scaling of a matrix X, we may use its condition number. If 
the matrix is poorly scaled, then its condition number is large. Minimizing the condition 
number therefore improves the scaling. For the above matrix, its condition number (in 
the Frobenius norm), defined as ||if p ||i7||(if p ) _1 || j p, is equal to 1158.2. If we choose p = 
y/l. 0000/66. 840 = 3.5000 • 10 _1 , the scaled characteristic polynomial has unit constant 
and highest coefficient, and the resulting scaled Hermite matrix reads 



SH P S 



163.48864 


151.37636 






-2445.0754 


163.00225 



151.37636 


2567.0923 




with 



S = diag(pVy,l). 
The Frobenius condition number of SH P S is equal to 32.096. 





163.00225 


150.92600 



Whereas this simple scaling strategy with one degree of freedom may prove useful for 
small-degree polynomials and small-size Hermite matrices, a more sophisticated approach 
is required for larger instances. 



4 Hermite matrix in Lagrange basis 

In this section we show how the Hermite matrix can be scaled by an appropriate choice of 
polynomial basis. Moreover, this basis allows for a straightforward entrywise construction 
of the Hermite matrix. 



4.1 Distinct interpolation points 

Consider n distinct interpolation points Uj G C, i = 1, . . . , n, and define the j-th Lagrange 
polynomial 

n 

U — Ui 



/,(«)= n 



, ,.Uj-Ui 



5 



which is such that lj(v>j 
1 



1 and lj(u. 



u 

2 



U 



u 



n— 1 



;) = if i 7^ j. In matrix form we can write 
1 • ■ ■ 1 



Ui 



U, 



u 



n-l n-l 



U 



n-l 





h 






h 






h 


» 




_ In 





(3) 



where V u is a Vandermonde matrix. Given a univariate polynomial q(s) with real coeffi- 
cients, define 

a(u) = Imq(jw) ,^ 
b(u) = Req(jw) 

as its imaginary and real parts on the imaginary axis, respectively. In the following, the 
star denotes transpose conjugation and the prime denotes differentiation, i.e. 

da{u) 



a!{u) 



du 



Theorem 4.1 When the interpolation points are distinct (i.e. m ^ Uj, i ^ j, i,j = 
1, . . . ,n), the Hermite matrix ofq(s) in Lagrange basis, denoted by H L , is given entrywise 
by 

aju^bjuj) - aju^bju*) 

11 U i / 



'It ■ 



UJ — Un 



a'(u*)b(uj) — a(uj)b'(u*) otherwise, 



for all i, j = 1, . . . , n. 

Proof Let us express the Bezoutian of a and b as a bivariate quadratic form 
a(u)b(v) — a{v)b{u) 



u — v 



1 


* 


1 


V 




u 


y n-l 




u n-l 



where H p is the Hermite matrix of q in the power basis. Recalling relation 
Bezoutian becomes 

a(u)b(v) — a(v)b(u) 



the 



l(v)*VJH i 'VJ{u) = l(v)*H L l(u) 



u — V 



so that the Hermite matrix of q in the Lagrange basis can be expressed as 

H L = V*H P V U . 

By evaluation at n distinct interpolation points Ui and v~, H L is given entrywise by 



K 3 



"(OKvj) - a{ Vj )b(u*) 



(5) 



Now let u* — > Vj for all i,j = l,...,n. After adding and subtracting a(vj)b(vj) to the 
numerator of (jSJ), we find 

H t L d = a'(u*)b( Uj ) - a( Uj )b'(u*), 
using a limiting argument. □ 
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4.2 Repeated interpolation points 

Let us define the bivariate polynomials 

d l+j ~ 2 f a(u)b(v) - a(v)b(u) 

Ci ' j[u,v} : ~ ^*=w=i v ^ 

for all i, j = 1, . . . , n and denote by 



{k) d a(u) 
W ~ du k 



the /c-th derivative of univariate polynomial a{u). 

Lemma 4.2 When the interpolation points are all equal (i.e. Ui = u,j for all i,j 
1, . . . ,n), the Hermite matrix of q(s) in Lagrange basis is given entrywise by 

(< - - 1)! *«»^> 
H ij ■= \ , 

, ^ wlse ' 

for all i,j = 1, . . . , n. 



Proof The proof of this result follows along the same lines as the proof of Theorem 14.11 
with additional notational difficulties due to higher-order differentations. □ 

Example 4.3 Let us choose n = 3 equal interpolation points (u\ = U2 = = x G M.). 
According to Lemma I4.2[ H L has the following entries: 

ttL a'(x)b(x)— a(x)b' (x) 

-"ll — li 

ttL (x)b(x)-a(x)b( 2 ) (x) 

n i2 — 2' 

ttL a^(x)b(x)-a(x)b( z ){x) 

""13 _ 3! 

(x)b'(x)— a'(x)b' 2 ' (x) . a' 3 ' a(x)b' 3 ' (x) 

-"22 ~~ 2> ' 3> 

rrL Q (3) (aQfe'^Q-a^x)^ 3 ) (a:) . agj (x)b(x)-q(x)fe( 4 ) (a) 

-"23 ~~ 3' + 4' 

rri q( 3 > (aQfe( 2 )(a)-q( 2 )(x)b( 3 ) (x) qW (^'(aQ-q'Qr)^ 4 ) (x) q< 5 > (x)b{x)-a(x)b^ (x) 

-"33 — 3!2! + 4! + 5! 

Based on Theorem 14.11 and Lemma 14.21 we leave it to the reader to derive entrywise 
expressions for the Lagrange basis Hermite matrix in the general case when only some 
interpolation points are repeated. 

In the remainder of the paper we will assume for notational simplicity that the interpo- 
lation points are all distinct. 
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4.3 Scaling 



Corollary 4.4 Let the interpolation points be (distinct) roots of either a(u) or b(u), as 
defined in Op. Then the Hermite matrix of q(s) in Lagrange basis is block diagonal, 
with 2x2 blocks corresponding to pairs of complex conjugate points and lxl blocks 
corresponding to real points. 



Proof From Theorem 14.11 all the off-diagonal entries of H L are given by 

a{u*)b{ Uj ) - aju^bju*) 

K - Uj 

when interpolation points Ui and Uj are not complex conjugate. Both terms a{u*)b{uj) 
and a(uj)b(u*) are equal to zero in since the interpolation points are the roots of either 
a{u) or b{u). The diagonal entries are a'{u*)b{uj) — a{uj)b'{u*) since it is assumed that 
interpolation points are distinct. Therefore this part of H L is 1 x 1 block-diagonal. 

When interpolation points and Uj are complex conjugate, there is only one non-zero 
entry which is equal to a'{u*)b{uj) — a{uj)b'{u*) and located in the off-diagonal 

entry, according to pairness. The diagonal entries of this case are equal to zero by virtue 
of equation ([6]). Therefore this part of H L is 2 x 2 block-diagonal. □ 

From Corollary 14.41 it follows that we can easily find a block-diagonal scaling matrix S 
such that the scaled Lagrange Hermite matrix 

H s = SH L S 



has smaller condition number. Nonzero entries of S are given by 

-i 



H 



whenever Hf- is a nonzero entry (i, j) of H L . 



Example 4.5 As an illustrative example, consider problem NN5 in [U] 
characteristic polynomial is 



The open-loop 



q(s) = s 7 + 10.171000s 6 + 96.515330s 5 + 458.42510s 4 

+2249.4849s 3 + 1.2196400s 2 - 448.72180s + 6.3000000. 



The Hermite matrix in power basis has the following entries: 



H 2,2 

H 5,5 
H 6fi 



-2826.9473 

608.04658 

-14719.034 

-4570.2494 

-4687.9634 

1026532.4 

21366.759 

523.23232 



H l,3 

H 3,7 

H 5,7 
H 7,7 



-14171.755 

-6.3000000 

206313.38 

209056.94 

1.2196400 

-22878.291 

-458.42510 

10.171000, 
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remaining nonzero entries being deduced by symmetry. Apparently, this matrix is ill- 
scaled. Choosing interpolation points itj as roots of a (it), the imaginary part of q(s) along 
the imaginary axis, we use Theorem 14.11 to build the Hermite matrix in Lagrange basis: 



diag (-2826.9473, 41032866 ■ 10 3 , 44286011 ■ 10 2 , 41032866 ■ 10 3 , 

22222.878 
22222.878 



44286011 • 10 2 , 



This matrix is still ill-scaled (with Frobenius condition number equal to 2.0983 ■ 10 7 ), but 
it is almost diagonal. Using an elementary diagonal scaling matrix S, we obtain 



H s = SH L S = diag (-1, 1, 1, 1, 1, 



1 

1 



which is a well-scaled representation of the Hermite matrix, with Frobenius condition 
number equal to 7. 



4.4 Target polynomial 

In our control application, let us introduce our main tuning tool which we call target 
polynomial, denoted by q(s). The target polynomial provides the interpolation points 
required to build well-scaled Hermite matrix in the SOF problem. These points are 
defined as in Corollary 14.41 as the roots of either the real or imaginary part of q(s) when 
evaluated along the imaginary axis. 

In the context of SOF design, the target polynomial may be either choosen as 

• a valid closed-loop characteristic polynomial ([I]) for a specific value of k, or 

• a polynomial with desired pole distribution for the closed-loop system. 

Furthermore, we invoke a continuity argument to observe that the condition and/or scaling 
of the Hermite matrix does not change abruptly in a neighborhood of a given target 
polynomial. 

Example 4.6 Consider again Example 12 .11 and let the target polynomial be an achievable 
closed- loop characteristic polynomial q(s, k) = det(sl — A — BKC), where 

K = [—4.3264 ■ 10 -1 , -1.6656, 1.2537 ■ HT 1 , 2.8772 ■ 10" 1 ] 

is a random feedback gain. The roots of the imaginary part of q(s, k) are chosen as 
interpolation points 

u = (0, ±60.847, ±16.007, ±9.2218, ±2.7034i) . 
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Here are two representative entries of the resulting Lagrange basis Hermite matrix: 

H^ 3 (k) = 9.4439251 ■ KT 1 + 1.9763715 ■ lCT 4 ^ - 8.9049916 ■ 10~ 4 A; 2 

- 8.6909277 ■ lO -3 /^ + 1.9212126 ■ lO^h 

+ 3.8300306 ■ lO" 9 *^ - 1.0276186 ■ lO -8 ^^ 
+ 3.3905595 • lO^hh - 3.4222179 • lO" 5 *^ 

- 3.8046300 • lO- 5 ^^ + 2.7420115 • 10" 9 £| 
+ 6.5442491 • lQ- 6 k 3 k 4 + 1.015195648 ■ lO^kj 



and 

Hl x {k) = -1.6918611*1 + 3.7288052 • lO" 1 ^ + 1.2286264 • Kr 2 JfeiJfe 3 . 

Comparing with the entries of the power basis Hermite matrix H p (k) given in Example 
12.11 we observe a significant improvement in terms of coefficient scaling. 



5 Numerical examples 



In this section, we present the benefits of Lagrange basis against power basis when solving 
SOF PMI problems found in the database COMPUb, see [9]. Even though Michal Kocvara 
and Michael Stingl informed us that an AMPL interface to PENNON is now available to 
solve PMI problems, in this paper for simplicity we consider only BMIs (i.e. quadratic 
PMIs) and the PENBMI solver (a particular instance of PENNON focusing on BMIs) 
under the YALMIP modeling interface, see [TU]. The numerical examples are processed 
with YALMIP R20070523 and PENBMI 2.1 under Matlab R2007a running on a Pentium 
D 3.4GHz system with 1GB ram. We set the PENBMI penalty parameter P0 by default 
to 0.001 (note that this is not the default YALMIP setting). 

As in [5J, the optimization problem to be solved is 

min fcjA — A 

s.t. H(k) y XI 

where H(k) is the Hermite matrix in power or Lagrange basis, fi > is a parameter and 
|.| is the Euclidean norm. Parameter \x allows to trade off between feasibility of the BMI 
and a moderate norm of the feedback gain, which is generally desirable in practice, to 
avoid large feedback signals. This adjustment is necessary in many examples. Indeed, 
the smallest values of ||fc|| are typically located at the boundary of the feasibility set, so 
the resulting closed-loop system is fragile and a small perturbation on system parameters 
may be destabilizing. 

PENBMI is a local optimization solver. Therefore, the choice of initial guess X.Q is 
critical. In most of the examples we choose the origin as the initial point. However this 
is not always an appropriate choice, as illustrated below. In addition to this, PENBMI 
does not directly handle complex numbers (unless the real and imaginary parts are split 
off, resulting in a real coefficient problem of double size), so we restrict the interpolation 
points to be real numbers. 
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As a result of the root interlacing property, the roots of real and imaginary parts of a 
stable polynomial are real (and interlacing). Owing to this fact, if we choose a stable 
target polynomial q(s) the resulting interpolation points are necessarily real. 

Example 5.1 Consider again problem AC4, with characteristic polynomial 

q(s, k) = s 4 + 150.92600s 3 + (130.03210 - 18.135000^ - 19612.500A; 2 )s 2 

- (1330.6306 + 19613.407A;i + 18322.789^)3 - (66.837750 + 980.62500£;i + 867.10818fc 2 ) 

and power basis Hermite matrix with entries 

H[ A = 88936.354 + 2615765.6A;i + 2378454.6fc 2 

+ 19233397A; 2 + 34974730^2 + 15887840A; 2 
H[ 3 = 10087.554 + 148001.81^ + 130869. 17A; 2 
H 2,2 = ~ 162937.14 - 2378239.7A;i + 23845311^ 

+ 355689. 13^ + 38500022 • lO 1 ^ + 35935569 • lO 1 ^ 
if£ 4 = 1330.6306 + 19613.407^ + 18322.789A; 2 
H^ 3 = 20955.855 + 16876.364^ - 2941713.4A; 2 
H[ A = 150.92600. 



Open-loop poles of the system are (2.5792 ,-5.0000 • 10~ 2 ,-3.4552, -150.00). If we 
define our target polynomial roots as (-5.0000 • 10~ 2 , -5.0000 • 10" 2 , -3.4552 , -150.00), 
keeping the stable open-loop poles and shifting the unstable open-loop pole to the left 
of the imaginary axis, our 4 interpolation points (roots of the real part of the target 
polynomial) are u = (±23.100 , ±4.9276 ■ 10~ 2 ) and the resulting Lagrange basis Hermite 
matrix has entries 

= 6.3432594 • 10" 1 + 3.1878941 • lO" 1 ^ - 17.462079fc 2 

+ 4.4822907A; 2 ± 4.4060140£;i£; 2 ± 4.1121739A; 2 
Hf 2 = - 3.7795293 • 10" 1 - 1.0581354*4 - 18.455273A; 2 

- 3.6574685 • 10~ 3 A; 2 - 4.4045142^2 - 4.1114926^ 
His = 2.4459288 + 36.285639A;i + 42.619220A; 2 

+ 7.8459729A: 2 + 189.06139^2 + 169.77324A; 2 
Hf A = 1.9481147 + 28.929191^ ± 12.037605fc 2 

+ 7.5224565A; 2 - 161.11487^^2 - 157.07808A; 2 
H$ 2 = 6.3432594 • 10" 1 + 3.1878941 • lO -1 ^ - 17.462079^ 

+ 4.4822907 • lO" 3 *; 2 + 4.4060140^2 + 4.1121739A; 2 
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H% 3 = 1.9481074 + 28.929083fci + 12.037568£; 2 

+ 7.5224134A; 2 - 161.11415^^ - 157.07737A; 2 
Hl A = 2.4459288 + 36.285639^ + 42.619220A; 2 

+ 7.8459729A; 2 + 189.06139^2 + 169.77324fc| 
H^ 3 = 659.47243 + 19434.408&1 + 18140.083fc 2 

+ 143181.92A; 2 + 267314.59^^2 + 124766. 18A; 2 
Hl A = 665.36241 + 19520.378^ + 17279.067fc 2 

+ 143169.06A; 2 + 253396. 76A;iA; 2 + 111775.41A; 2 
Hf A = 659.47243 + 19434.408^ + 18140.083A; 2 

+ 143181.92A; 2 + 267314.59^^2 + 124766.18A; 2 . 

Choosing the power basis representation with the orgin as initial point and trade-off 
parameter \x = 10 -5 , PENBMI stops by a linesearch failure and YALMIP displays a 
warning. However, we obtain a feasible solution A = 150.88 and K = [1.4181, —1.6809]. 
This computation requires 43 outer iterations, 433 inner iterations and 825 linesearch 
steps. On the other hand, in the Lagrange basis representation, the problem was solved 
with no error or warning, yielding A = 9.8287 ■ 10~\ K = [-5.0902 • 10" 2 , -2.0985 ■ 10~ 2 ] 
with 17 outer iterations, 100 inner iterations and 159 linesearch steps. 

We notice however that using the same trade-off parameter jj, for both representations 
is not fair since H p and H s have significantly different scalings. If we choose /i = 0.1 
for the power basis representation, no problem is detected during the process and we 
obtain A = 150.87, K = [8.0929 • 10" 2 , -1.6953 • 10" 1 ] after 26 outer iterations, 188 inner 
iterations and 238 linesearch steps. So it seems that the Lagrange basis representation 
becomes relevant mainly for high degree systems. This is confirmed by the experiments 
below. 

Consider the AC7, AC17, REA3, UWV, NN5, NN1 and HE1 SOF BMI problems of COMPLjb. In 
Table [T] we report comparative results for the power and Lagrange basis representations. 
As in Example I3.1[ the main strategy to choose the target polynomials (and hence the 
interpolation points) is to mirror the open-loop stable roots, and to shift the open-loop 
roots to, say —5.0000 • 10 _1 (any other small negative value may be suitable). We see 
that the behavior indicators of PENBMI are significantly better in the Lagrange basis, 
and the improvement is more dramatic for larger degree examples. More specifically: 

• for small degree systems like AC 17 there is only a minor improvement; 

• at the first attempt to solve the REA3 example strict feasibility was not achieved in 
the power basis, since A is almost zero. Therefore it was necessary to tune the \x 
parameter. Results of the second attempt show that the BMI problem was solved 
and the Lagrange basis computation was slightly less expensive than the power basis 
computation; 

• the underwater vehicle example UWV has two inputs and two outputs. However, 
because of cancellation of higher degree terms in the characteristic polynomial, the 
degree of the Hermite matrix is equal to 2 and we can use PENBMI on this problem; 
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• on open-loop stable systems such as UWV or AC17, the improvement brought by 
the Lagrange basis is less significant. Since the main purpose of our optimization 
problem is to minimize the norm of control gain, we observe that the Lagrangian 
basis is still slightly better than the power basis; 

• PENBMI is unable to reach a feasible point for examples NN5, NN1 and HE1, when we 
choose the origin as the initial point. Indeed, local optimization techniques seek an 
optimal point inside the feasible set in a neighborhood of the initial point. Therefore, 
achievement of the solver may be very sensitive to the initial point. When the initial 
point is defined heuristically or randomly, the improvement is significant for system 
NN5 in Lagrange basis. However, there is no improvement over NN1 and HE1, when 
we use this simple strategy to define the target polynomial. 



Table 1: PENBMI performance on SOF BMI problems 



system 


basis 




K 


out. 
iter. 


inn. 
iter. 


lin. 

steps 


K 


A 


AC7 

n = 9 


pow. 


1 


[0 0] 


27 


148 


167 


[1.1205 - 3.0946 • IO" 1 ] 


51.640 


Lag. 


IO" 5 


[0 0] 


15 


51 


67 


[5.7336 3.9995] 


3.6356 • 10 -i 


AC17 

n = 4 


pow. 


1 


[0 0] 


14 


65 


173 


[1.6619 • IO" 1 8.5782 • IO" 1 ] 


5.8306 


Lag. 


1 


[0 0] 


16 


36 


57 


-1.0855 • 10~* 1.5128 • 10 _1 


1.0459 


REA3 

n- 12 


pow. 


1 


[0 0] 


21 


28 


28 


- 1.0435 • IO" 5 - 2.2281 • 10" 4 


8.4187 • 10~ i3 


10" b 


[0 0] 


46 


458 


2460 


-43711 - 23491 


43787 


Lag. 


IO" 2 


[0 0] 


16 


48 


68 


- 4.2556 • 10 -i - 8.9973 • IO - * 


9.9105 • lO^ 1 


UWV 

n = 8 


pow. 


1 


[0 0; 0] 


13 


98 


188 


" -1.4319 • IO" 5 -2.6474 • 10 _e " 
-3.0817 • IO" 1 -5.6976 • 10" 2 


27.918 


Lag. 


1 


[0 0; 0] 


15 


65 


82 


-1.6755 • 10" i7 -6.9006 • 10" ia 
-3.6060 • IO" 8 -1.4851 • IO" 8 


1.0000 


NN5 

n = 7 


pow. 


1 


[10 5] 


29 


162 


300 


[12.382 9.0331] 


3.9610 • 10" 1 


Lag. 




[10 5] 


18 


45 


56 


[30.931 22.295] 


1.7652 • 10 _i 


NN1 

n = 3 


pow. 


i(H 


[0 30] 


15 


53 


59 


[7.9924 72.171] 


4.2238 


Lag. 


io~ 4 


[0 30] 


14 


49 


52 


[26.936 177.20] 


4.6019 


HE1 

n = 4 


pow. 


i 


[11] 


18 


73 


80 


[-1.5482 - 3.9063] 


34.359 


Lag. 


io- 1 


[11] 


18 


80 


87 


[-5.1376 11.589] 


32.168 



In Table [2] we show the influence of the target polynomial on the computational cost for 
the PAS system. Open-loop poles of the system are 

a = (0, 0, -9.5970 • 10 -1 , -36.646 ± 523.05i) 

and we choose alternative target polynomials with the following roots 

a x = (-5.0000 • 10~ 2 , -5.0000 • IO" 2 , -9.5970 • 10 -1 , -36.646 ± 523.05i) 
a 2 = (-1.0000 • IO" 3 , -1.0000 ■ 10~ 3 , -9.5970 ■ 10 _1 , -36.646 ± 523.05i) 
<7 3 = (0, -1.0000 • 10" 4 , -9.5970 • 10 _1 , -36.646 ± 523.05i) . 
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Table 2: Influence of target polynomial on PENBMI behavior 



system 


PAS degree 5 


basis 


power 


Lagrange 


Lagrange 


Lagrange 




1(T 3 


io- 8 


10~ 5 


io- 2 


roots 




o\ 






K 


[0 0] 


[0 0] 


[0 0] 


[0 0] 


out. iter. 


11 


19 


17 


15 


inn. iter. 


74 


27 


33 


29 


lin. steps 


194 


28 


44 


32 


K T 


-6.5390 • 10" 4 
-58.350 
-37.751 




-8.4106 • IO" 6 

-3.9048 
-9.9675- IO" 1 




-3.3369 • 10" 4 
-20.480 
-1.2157 




-8.6755 • 10~ 8 
-4.1040 • 10" 1 
-1.7471 • 10" 1 


A 


73.2917 


1.4901 • IO" 12 


8.1649 • IO" 3 


2.7241 • IO" 3 



One can easy to see that the computational cost is decreasing significantly when the point 
defining the target polynomial is getting closer to the PENBMI initial iterate. 

Consider the NN6 SOF BMI problem that was not solvable in the power basis, see [5]. 
Open-loop poles of the system are 

ct = (2.7303, 0, -7.2028 • 10~ 2 ± 60.804i, -1.0785 • 10" 1 ± 15.677i, -2.6764, -3.3000, -19.694) . 

The strategy to define the target polynomial is to change the unstable open-loop poles 
into slightly stable poles (shifting the real part to a small negative value). According to 
this strategy, our target polynomial has the following roots 

Gl = (-1.0000 -lQ- 3 ±i, -7.2028 ■ 10~ 2 ± 60.804i, -1.0785 ■ HT 1 ± 15.677i, 
-2.6764, -3.3000, -19.694) . 

The BMI SOF problem is solved with no error or warning in the Lagrange basis, yielding 
A = 8.8487- 10 _1 , K = [1.3682, 4.8816, 44.959, 59.016] with 17 outer iterations, 80 inner 
iterations and 138 linesearch steps, using the orgin as initial point and trade-off parameter 
fi = IO" 5 . 

6 Conclusion 

The Hermite matrix arising in the symmetric formulation of the polynomial stability 
criterion is typically ill-scaled when expressed in the standard power basis. As a conse- 
quence, a nonlinear semidefinite programming solver such as PENNON may experience 
convergence problems when applied on polynomial matrix inequalities (PMIs) coming 
from benchmark static output feedback (SOF) problems. In this paper we reformulated 
Hermite's SOF PMI in a Lagrange polynomial basis. We slightly extended the results of 
[T2] to use polynomial interpolation on possibly complex and repeated nodes to construct 
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the Hermite matrix, bypassing potential numerical issues connected with Vandermonde 
matrices. In our control application, a natural choice of Lagrange nodes are the roots of 
a target polynomial, the desired closed-loop characteristic polynomial. 

The idea of using the Lagrange polynomial basis to address numerical problems which 
are typically ill-scaled when formulated in the power basis has already proven successful 
in other contexts. For example, in [2] it was shown that roots of extremely ill-scaled 
polynomials (such as a degree 200 Wilkinson polynomial) can be found at machine pre- 
cision using eigenvalue computation of generalized companion matrices obtained by an 
iterative choice of Lagrange interpolation nodes. In [T3] the fast Fourier transform (a par- 
ticular interpolation technique) was used to perform spectral factorization of polynomials 
of degree up to one million. Another example of successful use of alternative bases and 
high-degree polynomial interpolation to address various problems of scientific computat- 
ing is the chebfun Matlab package, see pQ. Even though our computational results on 
SOF PMI problems are less dramatic, we believe that the use of alternative bases and 
interpolation can be instrumental to addressing various other control problems formulated 
in a polynomial setting. 

Appendix: Matlab implementation 

A Matlab implementation of the method described in this paper is available at 



Our implementation uses the Symbolic Math Toolbox and the YALMIP interface. It is not 
optimized for efficiency, and therefore it can be time-consuming already for medium-size 
examples. 

Let us use function hermitesof with its default tunings: 

» [A,B1,B,C1,C] = COMPleibCNNl') ; 
» A,B,C 
A = 

1 

1 
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homepages . laas . f r/henrion/sof tware/hermitesof . m 



B 





1 



C 



» [H,K] 
Quadrat i 




-1 



c 



5 -1 
-1 

: hermitesof (A, B,C) 

matrix variable 3x3 (symmetric, real, 2 variables) 
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Linear matrix variable 1x2 (full, real, 2 variables) 

Here are some sample entries of the resulting Hermite matrix 
» sdisplay(H(l,l)) 

-0 . 6168744435*K (2) -0 . 2372594014*K ( 1 ) *K (2) +0 . 04745188027*K (2) ~2 
» sdisplay(H(3,2)) 

0.3019687672*K(1)+0.01984184931*K(2)-0.009656748637*K(1)*K(2) 
+0 . 0003260141644*K(2) "2+0 . 04013338907*K(1) ~2 

For this example, the Hermite matrix is quadratic in feedback matrix K. This Hermite 
matrix is expressed in Lagrange basis, with Lagrange nodes chosen as the roots of the 
imaginary part of a random target polynomial, see the online help of function hermitesof 
for more information. In particular, it means that each call to hermitesof produces 
different coefficients. However these coefficients have comparable magnitudes: 

» [H,K]=hermitesof (A,B,C) ; 
» sdisplay(H(l,l)) 

-0 . 9592151361*K(2)-0 . 3689288985*K(1) *K(2)+0 . 0737857797*K(2) "2 
» sdisplay(H(3,2)) 

6 . 702455704*K(l)+0 . 5150092145*K(2)-0 . 3440108908*K(1) *K(2) 
+0 . 0184651235*K(2) "2+1 . 258426367*K(1) "2 

The output of function hermitesof is reproducible if the user provides the roots of the 
target polynomial: 

» opt = [] ; opt. roots = [-1 -2 -3]; 
» [H,K]=hermitesof (A,B,C,opt) ; 
» sdisplay(H(l,l)) 

-0.196969697*K(2)-0.07575757576*K(1)*K(2)+0.01515151515*K(2)~2 
» sdisplay(H(3,2)) 

0.2*K(1)-0.01818181818*K(2)-0.01212121212*K(1)*K(2) 
+0 . 0007575757576*K(2) "2+0 . 04166666667*K(1) "2 

The Hermite matrix can also be provided in the power basis: 

>> opt = [] ; opt. basis = 'p' ; 
» [H,K]=hermitesof (A,B,C,opt) ; 
» sdisplay(H(l,l)) 
-13*K(2)-5*K(1)*K(2)+K(2)~2 
» sdisplay(H(3,2)) 


For more complicated examples, the Hermite matrix H is not necessarily quadratic in K: 
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» [A,B1,B,C1,C] = COMPleib('NNl') ; 
>> size(B), size(C) 
ans = 

5 3 

ans = 

3 5 
» [H,K]=hermitesof (A,B,C) 

Polynomial matrix variable 5x5 (symmetric, real, 9 variables) 
Linear matrix variable 3x3 (full, real, 9 variables) 
>> degree (H) 
ans = 

5 
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